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1. EXECUTIVE SUMMARY 


The report concludes that It is possible to realize very substantial 
savings In the analysis part of the computational process In Integrated 
design and analysis of aircraft structures through the use of adaptive 
finite element technology based on the p~verslon of the finite element 
method. 

1.1 Background ; 

Adaptivity Is a procedui.e for efficient reduction of error on the 
basis of data already computed. The procedure Is comprised of two 
parts: determination of the sources of error In the computed data and 

estimation of their magnitude, and improvement of th'?. accuracy of approxi- 
mation by changing the number and/or distribution of the degrees of 
freedom. 

Currently there are only two adaptive finite element computer 
programs In existence: FEARS, developed at the University of Maryland 

and COMET, developed at Washington University. The two programs are 
based on completely different approaches: whereas FEARS Improves 

solution accuracy by selective mesh refinement, COMET upgrades the 
displacement approximation over finite elements through the addition of 
progressively higher order hierarchic shape functions.* 

This project involved some of the principal developers of both 
COMET and FEARS: Professor Ivo Babuska of the University of Maryland 

served as consultant to the project group at Washington University. 


*The term "p-verslon” Is defined In Section 3 (p. 9). 
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1.2 Main Results 

Adaptivity, based on the p-verslon of the finite element method 
will Impact on Integrated design and analysis In two ways: (a) substantial 

reduction In the cost of analysis and (b) Increased reliability of the 
computed data. These points are discussed In the following: 

1.2.1 Cost Reduction 

The most Important opportunity for cost reduction Is offered by the 
fact that the p-verslon exhibits faster rate of convergence than the 
conventional h-verslon based on uniform or quaslunlform meshes. Typi- 
cally, for comparable accuracy (say one to five percent error In strain 
energy) , the number of degrees of freedom required In the p-verslon Is 
only one-fifth to one-tenth the number required In the h-verslon. This 
Is especially Important In design optimization where a large number of 
analyses must be performed. The required computer time (cost) Is roughly 
proportional to the square of the number of degrees of freedom. Thus 
to reduction In the number of degrees of freedom results In to 

reduction In the required computer time In each analysis cycle. It 
Is noted that when the geometry Is complicated and a large number of 
elements Is required, the projected savings are similar due to the fact 
that the accuracy of analysis Is most strongly Influenced by the number 
and type of geometric constraint conditions (singularities) which tend 
to Increase with the geometrical complexity of the problem. 

Additional savings can be realized from optimizing the distribution 
of the degrees of freedom. This requires that the relative contribution 
of each element to the total error of analysis be known. The degrees of 
freedom are optimally distributed when the contribution of each finite 
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element to the total error of approximation Is the same. (This Is the 
basis for adaptive mesh refinement In FEARS, however the same Indicators 
cannot be used In the p-verslon) . Indicators were developed for the 
p-verslon for problems In two-dimensional elasticity under the present 
project. 

The savings that can be realized from optimal distribution of the 
degrees of freedom are highly problem dependent: In those cases where 

the number of singularities (reentrant corners, stiffener connections, 
changes In support conditions) Is large In relation to the number of 
finite elements, the savings will be small. In those cases where the 
number of singularities Is small In relation to the number of elements, 
the savings can be very substantial. 

Computational experiments were conducted with the objective to 
determine whether still further cost reductions can be realized through 
the use of the adaptive process In the following way: The number of 

degrees of freedom (hence the accuracy of approximation) was Increased 
gradually as the extremum was approached by the optimizer. The results 
of these experiments Indicate that the savings are only marginal. The 
reason Is that In statically Indeterminate structures the local optima 
may shift considerably as the accuracy of analysis Is Increased. Thus 
the optimizer may not be approaching the correct extremum value Initially, 
which negates the advantage gained from operating with fewer degrees of 
freedom at the beginning. 

A quantitative statement of error tolerance requirements generated 
In each optimization cycle by the optimizer could, conceivably, result 
In more substantial savings from this approach. 
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1.2.2 Reliability 

For obvious practical reasons, the finite element model should 
represent the prototype reliably. In other words, the accuracy of 
finite element analysis should not be sensitive to Input data, such as 
material properties, type of boundary conditions and loading. In order 
to verify that a finite element solution Is In fact correct, the user 
must ascertain that the solution Is in the asymptotic range either with 
respect to h -*■ 0 or p . -»■ <n . (h Is the maximum size of finite 

elements, is lowest polynomial order) . Because It can be 

extremely costly to verify that a given finite element solution Is 
Indeed In the asymptotic range with respect to h ->-0, such verifies- 
tlon Is almost Invariably omitted in practical analyses based on the 
h-verslon of the finite element method. In a number of cases, however, 
the accuracy of approximation shows extreme sensitivity to material 
properties and the kind of plate or shell theories used In the h-verslon, 
but not so In the p-verslon. Thus the p-verslon improves reliability In 
two ways: by making It simple to check whether a solution Is In the 

asymptotic range with respect to ^ by showing almost no 

sensitivity in the accuracy of analysis to variations In Input data. 

1.3 Recommenda t Ions 

In order to minimize computational costs associated with Integrated 
design and analysis processes, the following are recommended. 

(1) The analyzer should have p-verslon capabilities. It Is now a 
well established fact that the p-verslon of the finite element Is substan- 
tially more efficient than the conventional (h-verslon) In the three 
most costly areas of the analysis process: Input data preparation and 

modification, computer time, and verification of computed results. 
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(11) Because the cost of analysis very strongly depends on error 
tolerance n qulrements, the maxlmuin acceptable relative error should be 
determined by the optimizer for all constrained variables and communicated 
to the analyzer at each computational cycle. 

(Ill) Development of direct error estimation techniques for all 
constrained variables should be vigorously pursued. This Is a research 
problem rooted mostly In applied mathematics. The groundbreaking work 
has been completed, the existence of certain kinds of estimators has 
been proven. The objective of future research should be to extend the 
work for all constrained variables (deflections, moments, vibration 
frequencies) and problem types (stiffened plates and shells) encountered 
In structural synthesis. 

(Iv) Until direct error estimation techniques are available. 

Indirect methods should be used. This Involves the establishment of 
benchmark problems, such as those presented In this report, and develop- 
ment of correlations between the type of problem and the smoothness of 
approximating function on one hand, and relative error In the quantities 
of Interest on the other. 

(v) For structural synthesis Involving plates and shells finite 
element models based on the Relssner-Mlndlln theory are recommended. 

It has been demonstrated that the Relssner-Mlndlln theory Is capable 
of approximating the corresponding elasticity solutions to within one 
and three percent relative error In energy for a wide range of thick- 


nesses 
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INTRODUCTION 


There has been very substantial progress In the development o£ 
structural optimization methods In recent years. Unfortunately, their 
acceptance In design practice has been rather slow because efficient 
Integration of general purpose structural analysis systems with general 
purpose optimizers has not yet been achieved. An approach, which promises 
to remove this obstacle, was proposed In reference [1]. In this approach, 
a programming system consisting of a general purpose analysis program, a 
general purpose optimization program and two problem-dependent interface 
programs Is constructed. A specific implementation, called PROgramming 
^stem for Structural Synthesis (PROSSS) has been documented in [2]. 

In structural synthesis the structure Is modified and analyzed 
repeatedly until the optimality criteria are satisfied. The main cost 
source In the process Is the computer time required for structural 
analysis. Thus, to improve the efficiency of structural synthesis. It 
Is necessary to reduce the cost of analysis. 

Important developments have occurred In this area. In particular, 
the p-verslon of the finite element method has been developed and con- 
vergence rate theorems have been established for It. These theorems 
Indicate that In most problems of practical Interest the number of 
degrees of freedom required for a given level of precision Is substan- 
tially smaller In the p-version than In the conventional h-verslon which 
Is based on uniform or quasl-unlform mesh refinement. Also, it has been 
demonstrated In the course of the present project that the accuracy of 
analysis In the p-version Is not sensitive to input data such as material 
parameters and the types of plate and shell theories used to construct 
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the elemental stiffness matrices and load vectors. Techniques, such as 
reduced Integration, which Introduce uncertainties concerning the quality 
of approximation, are not needed In conjunction with the p-version even 
If Poisson's ratio Is as large as 0.4999 (near incompressibility) or 
when shear deformation Is accounted for In the analysis of thin plates. 

The convergence rate theorems are a priori estimates of error l.e. 
the estimators do not require that the computations be actually performed. 
These theorems are useful for making general comparisons between alternative 
modeling strategies. It Is desirable to obtain more accurate error 
estimates however, applicable to specific approximations. Error estimators 
based on computed results are called a posteriori estimators. It Is 
reasonably certain that reliable a posteriori error estimators can be 
developed for the p-verslon of the finite element method, but some 
theoretical problems remain. The Importance of such estimators cannot 
be overstated. They are essential for attaining the goal of computations 
which Is either to achieve some desired level of precision at minimum 
overall cost or to obtain the best possible approximation for a fixed 
cost with a reliable estimate of the error of approximation. The cost 
cannot be minimized If, due to lack of a reliable error estimator, the 
analyst Is forced to employ more degrees of freedom than necessary in 
order to be "on the safe side". It Is self evident that the quality of 
an approximation cannot be accurately assessed if means for error estima- 
tion are not available. For these reasons, a great deal of attention 
was given to a posteriori error estimation in the present project. 
Significant progress has been made in establishing that reliable local a 
posteriori error Indicators exist In the p-verslon at least for one and 
two-dimensional elliptic boundary value problems. 



In this report we first sunanarize the main mathematical theorems 
which establish the p-version of the finite element method. The theorems 
give a priori error estimates. We then discuss practical considerations 
resulting from the theorems with respect to computational efficiency. 

The discussion is supported by examples which are suitable for benchmark 
comparisons with other computer codes. We also present examples of 
structural optimization. Finally, the development of a posteriori error 
Indicators is discussed. 
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3. DEFINITIONS AND PRELIMINARIES 

There are two basic convergence processes in finite element analysis. 
The conventional finite element convergence process in which the number 
of Interpolating (shape) functions are fixed for each element and the 
mesh Is refined In such a way that the maximum diameter of elements, h, 
approaches rero Is called h-conver,tence . In the other basic convergence 
process, called p-convergence . the number and distribution of finite 
elements Is fixed and the interpolating functions, which are complete 
polynomials of order p, optionally supplemented with other types of Inter- 
polating functions, is progressively Increased. 

There are very substantial differences between the two convergence 
processes in numerical performance and computer Implementation. For this 
reason we refer to the conventional computer implementation of the finite 
element method in which Improvement of accuracy Is achieved by mesh 
refinement as the h-verslon of the finite element method. Computer 
implementations In which Improvement of accuracy Is achieved by Increasing 
the polynomial order over a fixed mesh are referred to as the p-verslon 
of the finite element method. 

In the Interest of programming and computational economy, the p-version 
utilizes exactly and minimally confovmlng hierarchic finite elements . 
Hierarchic finite elements have the property that the set of basis 
functions of an element of order p is a subset of the basis functions of 
all higher order elements of the same kind. Consequently, the stiffness 
matrix of an element Is embedded in the stiffness matrices of all higher 
order elements of the same kind. 

Exact and minimal conformity Is important for the following reasons: 
Finite element approximation is based on minimizing the total potential 
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energy 7t(u) which, In the absence of body forces. Initial stresses and 
strains that would not affect the following discussion, can be written 
as 


ir(u) 


1 

2 



([D]u)^[E]([Dlu)dV - 



T 

u 


T ds 


(3.1) 


in which [D] is the differential operator matrix that defines the strain- 
displacement (or equivalent) relations, u is the displacement vector 
field, [E] is the matrix that defines the stress-strain relations 
(Hooke's law or equivalent), T is the surface traction vector. The 
first integration is over the entire volume of the structure, the second 
is over that part of the surface for which surface tractions are specified. 

The theorem of minimum potential energy is usually stated in the 
following way: "Of all displacements satisfying the given boundary 

conditions, those which satisfy the equilibrium equations make the 
potential energy an absolute minimum". (See, for example, reference 3). 

In finite element analysis we are Interested in the converse of this 
theorem, namely, we wish to minimize the potential energy functional 
tt(u) over some set of admissible trial functions which will then 
approximately satisfy the equilibrium equations. It is well known that 
the admissible functions must satisfy the principal boundary conditions 
[3]. It is almost self evident that the set of admissible functions 
should be as large as possible, in order to permit us to make the 
error | | (measured in some suitable norm, to be discussed later) 

as small as possible. To explain: if we arbitrarily restricted the 

class of admissible trial functions, we certainly would not have 
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Improved our chances for minimizing the error In Che most efficient way 
possible. For this reason It Is desirable that the Interelement con- 
tinuity requirements be minimally satisfied: Any continuity requirement 

over the minimum Imposes an unnecessary restriction on the class of 
admissible trial functions. 

Let us now examine the consequences of not satisfying the minimal 
continuity requirements exactly. The continuity requirements arise from 
the criterion that admissible trial functions. In our case must 

result In finite strain energy. Thus, any trial function for idilch the 
first Integral In eq. (3.1), called the strain energy Is Infinitely 

large, la not admissible. Therefore Is admissible If It satisfies 
the principal boundary conditions and Its first partial derivatives 
are square Integrable on the solution domain. In mathematical notation, 

Is admissible If e The notation reminds one of the 

essential characteristics of Each component of vector belongs 

to the space of functions H, which are defined on the solution domain Q, 
satisfy the principal boundary conditions on the boundaries of n, 
denoted by F, and have square Integrable first derivatives. 

It Is not difficult to show that If the trial function Is not 
continuous from one finite element to the next then U(^^) Is unbounded. 

For example, let us consider the Interelement boundary between element A 
and element B In Fig. 3.1. 

For the sake of simplicity let us assume that Ug(o,y) - u^(o,y) > 0 
in the range y^^ < y < y 2 * Let M^(e) represent the least upper bound 
(supremum) of u^ in the region - e/2 < * < 0» y^^ < y < y 2 Let M^Ce) 
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represent the greatest lower bound (Infimum) of In the region 
0 < X < e/2, 7]^ 1 y < ^ 2 * ^ ^ ® arbitrarily. Then: 


11m 

e-»-0 



e/2 




2 


dxdy = 


(3.2) 



Fig. 3.1 

Notation for equation 3.2 



-13- 


Because ■ ) occurs in the strain energy expression, the 
strain energy Is not defined for discontinuous u^^. Nonconforming 
elements are not admissible. This does not mean that nonconforming 
elements are not capable of yielding "good results". It does mean that 
the numerical performance of nonconforming elements is problem»dependent : 
It Is possible to load nonconforming elements in such a way that they 
deform without absorbing strain energy. Consequently the question of 
whether the results obtained by means of nonconforming elements are 
reliable In any given situation must be addressed. We avoid this problem 
by employing exactly and minimally conforming finite elements and the 
displacement formulation. 

References 4 to 8 present detailed descriptions of exactly and 
minimally conforming hierarchic finite elements. 

3.1 Rates of convergence; uniform or quasi-uniform mesh refinement 

and p-dlstrlbutlon 

We shall now summarize the main results concerning the rates of 
convergence of the finite element method. For further details, references 
9, 10, 11 should be consulted. 

The relationship between error and the number of degrees of 
freedom In finite element analysis Is governed by a property of the 
approximated function, called "smoothness". With some restrictions to 
be noted later. It Is generally true that the smoother the approximated 
function, the fewer the number of variables needed to achieve a given 
level of precision. Angular corners at external boundaries tend to 
reduce smoothness. In twu-dlmenslonal elasticity, for example, the 
approximated displacement vector function Is of the form: 
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u - r“ ^(0) + G(r,0) (3.1.1) 

In which r and 0 are polar coordinates with the origin at the corner; 
a > 1/4 Is a smoothness parameter. (The smoothness of u Increases with 
a). Its value depends on the comer angle axid the boundary conditions 
Imposed at the comer; £ and G are smoother functions than u In the 
neighborhood of r > 0. The values of a for plane elastic and plate 
problems were given by Williams [12]. 

When the number of degrees of freedom, N, Is Increased through 
uniform or quasl-unlform mesh refinement (h-verslon) , the strain energy 
of the error Is bounded by: 

In which U Is the strain energy; c Is a constant which depends on the 
order p of the polynomial approximation, the value of elastic constants 
and the geometry of the mesh. 

A sequence of mesh refinements Is quasl-unlform If In the refinement 
process the ratio of largest element diameter to the smallest tends to 
a finite number. 

When the number of degrees of freedom Is Increased through Increasing 
the order of polynomial approximation (p-verslon), the strain energy of 
the error Is bounded by: 


(3.1.3) 
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In which k is a constant which depends on the elastic constants, the 
geometry of the mesh and the distribution of polynomial orders over the 
mesh. This relationship was established only very recently [10] . 

The exponents of N In eqs. 3.1.2 and 3.1.3 are called the asymptotic 
rates of convergence . When log U(^-Uj,g) Is plotted against log N, then, 
according to eqs. 3.1.2 and 3.1.3, for sufficiently large N, a straight 
line Is obtained whose slope Is the asymptotic rate of convergence. When 
the energy of the error of a given finite element solution Ic plotted 
against N on a log-log scale, and It falls on the straight line portion, 
then the finite element solution Is In the asymptotic range . 

Naturally, we want our finite element solution to be reliable. 

This means that the error should not depend significantly on Input 
parameters such as the choice of mesh, polynomial order or material 
properties but rather It should depend only on the function being 
approximated. For this reason a finite element solution should be In 
the asymptotic range and when the p-verslon or the h-verslon with 
uniform or quasl-unlform mesh refinement Is used. It should be governed 
by the parameter a. 

The foregoing summary touches on some of the most Important charac- 
teristics of the finite element method. In order to underline and 
Illustrate the main points, a simple example Is presented. 

3.2 Example; Parabollcally Loaded Square Panel 

The parabollcally loaded square panel Is Illustrated In Fig. 3.2.1. 
The approximated displacement vector field u Is exceptionally smooth 
In this case. Specifically, the value of a (defined In eq. 3.1.1) Is 
2.542, whereas In practlal problems a Is usually between 0.5 and 0.75. 




Fig. 3.2.1 

Farabollcally loaded square panel 
Typical finite element mesh 

The problem was chosen to demonstrate that in smooth problems the 
polynomial order controls the rate of h-convergence until p > a. 

The results for Poisson's ratio of 0.3 are shown In Fig. 3.2.2. 

For p * 1 and p ■ 2 the asymptotic range Is entered at about SO 
degrees of freedom In the h-verslon. The slopes for N > 50 (l.e. the 
rates of convergence) are governed by p, as predicted by eq. (3.1.2). 

For p * 3 the rate of convergencf s governed by a, again as predicted 
by eq. (3.1.2). Because of the exceptional smoothness of the approximated 
function, the relative errors are exceptionally small. The faster rates 
of p-convergence for the different meshes are shown by the dashed lines. 
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NUMBER OF DEGREES OF FREEDOM N 



In N 

Fig. 3.2.2 

Parabollcally loaded sqiiare panel 
Relative error in strain energy vs, N. Poisson's ratio: 0.3 

The results for Poisson's ratio of 0.4999 (near incompressibility) 
are shown in Fig. 3.2.3. For p « 1 we observe that the pre-asymptotic 
rate of convergence is extremely slow. In fact the asymptotic range 
is not entered until the roundoff limitations of most digital computers 
are reached. This example shows that when the solution is substantially 
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unchanged when the number of degrees of freedom Is Increased, does not 
necessarily mean that the solution "has converged". In this case It has 
not even entered the asymptotic range. 


NUMBER OP DEGREES OF FREEDOM (N) 
10 20 50 100 200 500 



2 3 4 5 6 7 

InN 


Fig. 3.2.3 

Farabollcally loaded square panel 
Relative error In strain energy vs. N. Poisson's ratio: 0.4999. 

For p * 2 the pre-asymptotlc rate of convergence Is slower than 
the asymptotic rate. For p * 3 the pre-asymptotlc and asymptotic rates 
are essentially the same. 
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For the p-verslon the much faster rate of convergence for the 
different meshes Is Indicated by the dashed lines. 

The poor convergence characteristics of the h-verslon in the case 
of nearly incompressible solids are well known (see, for example [13]). 
The reduced integration technique, which raises new questions relating 
to quality control, has been proposed. The p-verslon eliminates the 
need for reduced integration. 

The behavior of h-convergent approximations is similar in the case 

of thin plates and shells based on Mlndlin's theory: When shear is 

included in the strain energy expression, its significance in relation 

to the flexural strain energy progressively diminishes as the thickness 

is reduced. Mathematically, a penalty multiplier is applied to the 

shear strain energy which Increases, relative to the flexural strain 

_2 

energy, in proportion to t , t being the plate thickness. This is 
analogous to the case of nearly incompressible solids, in which the 
penalty multiplier is a function of Poisson's ratio and is applied to 
that part of the strain energy which is associated with volumetric 
strain. 

3.3 Rates of Convergence: Non-quasiuniform Mesh 

Substantial research effort has been devoted to finding ways to 
Increase the rate of convergence of the h-verslon of the finite element 
method. The most Important result of this effort is that it has been 
proven in [14] that there exist sequences of non-quaslunlform meshes 
such that the rate of convergence in the h-version is governed by the 
polynomial order p, not by the corner singularities [l.e. parameter a 
in Eq. (3.1.1)]. Thus with proper mesh refinement, the rate of conver- 
gence is: 
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h (3.3.1) 

rather chan Che race given by eq. (3.1.2) which Is valid for uniform and 
quaslunlform mesh refinements. The proper mesh refinement Is problem 
dependent, however. This means chat the meshes must be determined 
adaptively. The adaptive mesh refinement procedure Is based on error 
estimators. The estimators measure the error contribution of each 
element to the total error of approximation. The mesh Is then refined 
so as to make the error contribution from each element as nearly uniform 
as possible. The procedure has been Implemented In a computer code 
called FEARS (F^lnlce Element ^aptlve Research ^olver) at the University 
of Maryland [15]. 

The finding chat the asymptotic rate of h-convergence can be made 
Independent from the corner singularities and In fact arbitrarily large, 
since the polynomial order p can be chosen arbitrarily. Is certainly one 
of the most Important theoretical results concerning the convergence 
characteristics of the finite element method. It Is doubtful, however, 
that this result can be exploited for purposes of quality control In 
general purpose finite element computer programs. The algorithmic 
structure of adaptive finite element codes based on the h-verslon Is 
very complicated due to the fact chat the mesh changes at each adaptive 
step. The cost of data management operations can substantially exceed 
the cost of solution of the linear systems of equations. Furthermore, 
adaptivity In Itself does not alleviate the sensitivity of the h-verslon 
to Input parameters such as Poisson's ratio and plate thickness. 
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In the p-verslon the asymptotic rate of convergence cannot be 
made Independent from the corner singularities by adaptively 
selected p-distributions but it is possible to reduce the number of 
degrees of freedom without reducing the accuracy by optimizing the 
p-distribution. The basis for obtaining optimal p-distribution is 
similar to that used in the h-version: Error indicators are computed 

which measure the relative contribution of each element to the total 
error of approximation. The p-distribution is optimal when the error 
contribution of each element is the same. At present the relative 
error contributions can be e.^tlmated for plane elastic problems from 
the residuals obtained when the finite element solution is substituted 
into the Navler equations. This development will be discussed in 
detail in section 6. Error estimators, similar to those in [14] are 
not yet available for the p-verslon. 
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4. THE COST OF ACCURACY 

The goal of conputatlona in engineering practice In general and 
structural synthesis In particular Is either to achieve some desired 
level of precision at mlnlnun overall cost or to obtain the best 
possible approxlioatlon for a fixed cost with a reliable estloate of 
the ?*cror of approxlnation. 

The desired level of precision Is typically around 1 percent in 
strain energy and 1 to 3 percent In stress resultants (moments and 
shear forces); stress Intensity factors and vibration frequencies in 
linear analyses. It Is not meaningful to measure the error In terms 
of stress maxima because due to the almost Invariably present corner 
singularities the maximum stress computed on the basis linear elas- 
ticity Is Infinitely large. 

The overall cost of computations Is comprised of three major 
parts: Data preparation, computer time and verification of results. 

Because data preparation usually Involves large amounts of human 
effort. It Is generally regarded as the largest of the three cost 
Items. The cost of computer time Is diminishing In relation to the 
others, nevertheless in certain types of analyses, such as optimiza- 
tion and nonlinear analyses, computer time remains the controlling 
factor. State of the art (h-vcrslon) computer codes do not provide 
methods for verification of computed data which are accepted In 
general on the basis of the analysts' professional Judgement 
developed through experimentation with benchmark problems. Such 
Indirect methods of verification are not always reliable In the 
h-verslon, however. For example, we have seen in the case of nearly 
incompressible solids that no change In the strain energy with respect 
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to increasing the degrees of freedom by mesh reflnemenc did not guarantee 
the solution error to be small. The reason for this Is that the pre- 
asymptotlc rate of convergence of the h-version Is sensitive to Input 
parameters. In this case to Poisson’s ratio. 

We shall now consider how the costs of structural analyses may be 
reduced and the reliability Increased by means of existing or proven 
technology. 


4.1 Comparison of the h and p versions In the asymptotic range 

Let us first assume that finite element solutions were obtained 
by both the h- and p-verslons for a given problem and the approxima- 
tions are In the asymptotic range. Let us compare, under these condi- 
tions, the computer costs of Increasing accuracy In the computed 
strain energy by one significant digit, using uniform mesh refinement 
and p-dlstrlbutlon, and the costs Involved in structural synthesis. 

We note that the asymptotic convergence rate for the relative 
error In strain energy, e, defined as: 


iHu-ajj) 


(4.1.1) 


Is the same as asymptotic rates given by eqs. 3.1.2 and 3.1.3. Denoting 
the relative error corresponding to d significant digits of precision 
by e^, we have, by definition: 

«d 


(4.1.2) 
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The major variable cost in the solution of finite element problems 
Is associated with solving the corresponding system of simultaneous 
equations. In general: 

(4.1.3) 

in which C Is the cost of computer time and g Is a number between 1 

and 3. 3 depends on the bandedness or the front width of the linear 

system of equations and consequently on whether the h- or p-verslon Is 

used. Dependence on the mode of convergence Is quite weak, however, 

and can be ignored for the purposes of the present discussion. In 

two-dimensional applications the usual value of 3 Is close to 2. For 

example, the half band width B of the structure stiffness matrix 

corresponding to a uniform mesh on a square domain Is proportional to 
1/2 

N . The number of operations required for solving the linear system 

2 

of equations Is proportional to NB , hence 3*2. In computational 
experiments Involving p-convergent approximations 3 was found to be 
close to 2 also when Irons' frontal solution technique was utilized. 
Thus, although the stiffness matrix tends to be more fully populated 
in the p-verslon than In the h-verslon, sparse matrix solution tech- 
niques provide substantial reduction In the number of operations as 
compared with solvers that do not account for sparsity (3 * 3} . 

We remark that only the number of Internal modes Increases In propor- 
2 

tlon with p . The Internal modes are eliminated at the element level. 
The number of external modes, which represent connectivity among 
finite elements and must therefore be eliminated at the global level. 


Increases only in proportion to p. On the other hand, both the 
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-2 

external and Internal modes Increase in proportion to h . This tends 
to compensate for the more heavily populated stiffness matrices in the 
p-version. 

From the rates of convergence and cost formulas (eqs. 3.1.2« 3.1.3 
and 4.1.3) we find that the cost Is related to the error as: 

C ~ g-6/®ln(p,o) h-version (4.1.4) 

and 

C in the p-verslon (4.1.5) 

Denoting the cost corresponding to d significant digits of precision 
by C^, we obtain 

= j^Q0/®in(P»oi) h-version (4.1.6) 

^d 

and 

Q 

■ 10®^^“ in the p-version (4.1.7) 

In linear elastic fracture mechanics, for example, a ■■ 1/2. 

Assuming that 8 * 2, we find that the cost increases by the factor of 
10,000 for each additional significant digit in the h-verslon and by 


a factor of 100 in the p-version 
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For structural synthesis, the difference In computer costs between 
the h- and p-verslons based on uniform mesh refinement and p-dlstrlbutlon, 
can be estimated as follows: Let us consider a tjrplcal error versus 

degrees of freedom diagram, as shown In Fig. 4.1.1. The diagram 
Indicates that two extensions of an Initial solution (representing the 
minimum number of finite elements needed to define the geometry and 
the lowest polynomial order) were obtained, one by the h-verslon, the 
other by the p-verslon, until the relatve error e was reduced to some 


Log N 



Log e 


Fig. 4.1.1 

Typical relative error vs. number of degrees of freedom 
diagram representing refinement of an Initial solution by 
the h-verslon (assuming uniform mesh refinement) and the 
p-verslon (assuming uniform p-dlstrlbutlon) In the 
asymptotic range. 
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accep table level e . The corresponding degrees of freedom are denoted, 

respectively, by and N^. As long as the mesh refinement Is uniform 

or quaslunlform and singularities occur only at the vertices of finite 

elements, N, < N . Since the cost Is approximately proportional to 
h p 

the square of the number of degrees of freedom, the cost differential 
for large N, denoted by AC, can be estimated as: 

AC'wn(N^-Np) (4.1.8) 

In which n Is the number of times the analyses Is called In the synthesis 
process. We give specific examples of cost savings In structural 
synthesis In this report. We emphasize that the estimate represented 
by (4.1.8) Is valid for large N values only. 

4.2 Practical Considerations 

The foregoing analysis was based on the assumption that both the 
h- and p-verslons are In the asymptotic range. It has been observed 
by several Investigators, however, that In practical analyses the 
as}rmptotlc range Is very seldom entered In the h-verslon. This means 
that the solution error Is governed by Input parameters rather than 
the true solution of the problem. The pre-asymptotlc rate of conver- 
gence may be faster or slower than the asymptotic rate. In the 
p-verslon on the other hand the asymptotic range Is entered at low 
values of p and, very importantly, the point of entry Is substantially 
Independent from the Input parameters. 

In view of these facts, let us now ask the following question: 

Given that an analysis has been performed with a minimum number of 
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flnlte elements and low polynomial order, and it Is desired to Improve 
the accuracy of the solution; which Is the better strategy to follow: 
mesh refinement or Increasing the polynomial order? The answer must be 
based on two considerations: cost and reliability. Both considerations 

suggest that Increasing the polynomial order Is the better strategy. 

The costs of data preparation are certainly less In the case of the 
p-verslon since mesh refinement Is not Involved. The cost of computer 
time Is also less because fewer degrees of freedom are needed to achieve 
the same level of accuracy In the p-verslon than In the h-verslon. And, 
finally, the p-verslon Is more reliable than the h-verslon because the 
entry point Into the asymptotic range Is not significantly affected by 
the Input parameters. 

Compared on the basis of computer costs, the p-verslon has a very 
strong advantage over the h-verslon within the range of accuracies 
required In engineering analyses. The writers' experience suggests 
that even when very severe corner singularities are present, such as 
In linear elastic fracture mechanics, the relative error In strain 
energy can be reduced to under 1 percent by using coarse meshes and 
polynomial orders ranging from 6 to 8. For less severe comer 
singularities the polynomial orders are typically In the range of 3 
to 5 for 1 percent error. If, for any reason, the accuracy of analysis 
must be Increased such that the relative error Is substantially below 
1 percent, the best strategy Is to combine the h- and p-verslons by 
using strongly graded fixed meshes at comer singularities and In- 
creasing p until the desired level of precision Is reached. An 
example Is presented In [9]. Theoretically It Is possible to Increase 
both the accuracy and the rate of convergence beyond any limit by 
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optlmal mesh grading coupled with optimal p-dlstrlbutlons. The need 
for such techniques In practical applications Is not foreseen, however. 

The cost of verification of computed data In the p-verslon by 
means of repeated analyses Is further reduced by the hierarchic pro- 
perty of the shape functions: Suppose that an analysis Is to be 

repeated using selectively or uniformly Increased pol 3 momlal orders. 

Because the already triangulated stiffness matrix Is embedded In the 
new stiffness matrix, the elimination process must be continued only 
for the new rows and columns to complete the new analysis. This 
feature minimizes the marginal cost of verification by repeated 
analyses . 

4.3 Example: Double-edge cracked square panel 

We shall demonstrate on the basis of the double edge cracked 
square panel shown In Fig. 4.3.1 that engineering accuracy can be 
achieved with the p-verslon even when severe singularities are present. 

As In the case of the parabollcally loaded square panel, the problem 
has been analyzed for two different Poisson's ratios (v > 0.3 and 
V ” 0.4999) . The problem was also analyzed by the h-verslon. The 
results for v = 0.3 are shown In Fig. 4.3.2. It Is seen that a 
coarse mesh, consisting of only 8 elements, and p ranging from 3 to 6 
Is adequate for obtaining an approximation for which the relative error 
In energy Is between 1 and 5 percent. We note that the error In strain 
energy Is the logical measure for error because the stress Intensity 
factor Is related to the strain energy release rate. In fact the 
rate of convergence of the strain energy release rate has been demon- 
strated to be the same as the rate of convergence of the strain energy [16]. 
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The h-verslon, based uniform mesh refinement, requires a very 
large number of elements to reduce the relative error to 5 percent. 
When adaptive (non-quaslunlform) mesh refinement Is used (Fig. 4.3.3), 
the rate of convergence of the h-verslon Is governed by eq. 3.3.1. In 
this case p«l and 2o>l hence the same convergence rate Is predicted 
as for the p-verslon by eq. 3.1.3. 

When Poisson's ratio approaches 0.5, the point of entry Into the 
asymptotic range shifts toward very high N values for the h-verslon 
but not for the p-verslon. This Is shown In Fig. 4.3.4. The results 
are similar to those obtained for the parabollcally loaded square 
panel (Fig. 3.2.3). This degradation occurs In connection with the 
h-verslon even when adaptive mesh refinement Is used. 



Fig. 4.3.1 


Double edge cracked square panel 
Initial finite element mesh 


n)n - (n)n 



2 3 4 5 6 7 

In N 

Fig. 4.3.2 


Double edge cracked square panel 
Relative error in strain energy vs. N. 
Poisson's ratio : 0.3. 


RELATIVE ERROR (PERCENT) 








RELATIVE ERROR (PERCENT) 
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5. DIMENSIOHAL REDUCTION: VALIDITY OF PLATE BENDING MODELS 

A great deal of work has been devoted in the past twenty years to 
the development of plate bending finite elements. The papers on this 
subject are too numerous to mention but the problem itself and possi- 
ble solutions of the problem are adequately described in references 
17, 18, 19. Most of the work in this area has been concerned with 
the enforcement of continuity which is required for modeling 
plate behavior on the basis of Kirchhoff's theory. 

We approach the problem from a different point of view: Empha- 
sizing that the goal of computations is to approximate the solution 
of the equations of three dimensional elasticity over the plate 
domain, the thickness of which is usually much smaller than its 
other dimensions, we demonstrate that this goal is best achieved 
through a formulation which is based on the Relssner-Mindlln plate 
theory [20,21] and the p-version of the finite element method. 

This formulation has very Important advantages: 

(1) It is unnecessary to distinguish among "thin", "moderately 
thick" and "thick" plates because the numerical performance of the 
formulation is insensitive to plate thickness. (Property of robust- 
ness) . 

(11) The treatment of curved boundaries does not pose unusual 
problems because only C” continuity is required. C* continuity is 
preserved in transformations to curvilinear coordinates (such as 
in isoparametric transformation). 

(ill) The rate of p-convergence is very rapid. The number of 
degrees of freedom required for achieving engineering accuracy 
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(1 to 3 percent error in energy) Is generally small even when severe 
corner singularities are present. 

We begin our analysis with a review of the Relssner-Mlndlin plate 
theory and its finite element implementation. 

5.1 The Reissner-Mindlin Plate Bending Element ; 

For the sake of simplicity, let us consider an infinitely long 
plate strip of width L and thickness t loaded and supported in such 
a way that the problem can be viewed as a plane strain as well as a 
plate bending problem (fig S.1.1). Let us compare, under these 
conditions, the formulations of the Reissner-Mindlin plate theory 
and the three-dimensional theory of elasticity. 

The assumed displacement functions in the Reissner-Mindlin theory 

are: 

Uj^ - 

U2 - 0 (5.1.1) 

“3 “ '^^* 1 ^ 

The assumed displacement functions in the three-dimensional theory 
of elasticity are: 

00 

2i-l . . 

“1 “ 2. "*3 ♦2i-l^*l^ 

i-1 

U2 ■ 0 

00 

“3 ■ 2 

i-0 


(5.1.2) 
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In which the functions u^(x^,X 2 )i u^(x 2 iX 2 ) have been written as sums of 
products of powers of the transverse variable x^ and the field functions 
<^( 21 - 1 ) (x^) and <l'^21)^^l^* that the Relssner-Mlndlln theory 

retains only the leading terms of the series which spans the solution 
space, l.e. is capable of approximating the elasticity solution with 
arbitrary precision. 

The strain-displacement relations are the same In both cases: 

'ij ■ 2 <“l.j + “j.l^ 

The stress-strain relations differ significantly, however. In 
the case of elasticity: 

“ij ’ “ 'ij 

In which X and G are the Lame parameters. 

In the Relssner-Mlndlln plate theory the additional assumption 
Is Introduced that o^ • 0. We note that this assumption contradicts 
the assumption implied in eg. (5.5.1), namely that the transverse 
strain Is zero. Consequently, the stress-strain relationship is: 

°11 “ T~1 '^ll 
1-v 

vE 

®22 " , 2 ^22 
1-v 


(5.1.5) 
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0 


^13 


K£ 

£ 

2(l+v) 


In which R is the shear factor » to be discussed In Section 5*2. 

In both cases the strain energy expression is based on the fundamental 
expression: 


U 





(3.1.6) 


However, because of the differences In the stress-strain relationship, the 
strain energy expression of the Relssner-Mlndlln plate theory will not be 
the same as the strain energy expression of the theory of elasticity 
truncated so that only the first beams are retained for u^ and u^. 

The question arises whether it is reasonable to expect solutions of 
plate problems obtained by means of the Relssner-Mlndlln theory to remain 
close to corresponding solutions of the three dimensional theory of 
elasticity when the plate thickness Is large with respect to the other 
dimensions of the plate; with ret pect to the minimal wave length of 
oscillatory loadlxig or, more generally, with respect to the minimum 
wave length of the approximated function. This question has been Inves- 
tigated on the basis of examples presented In sections 5.3, 5.4 and 5.5. 
The conclusion Is that within the range of engineering accuracies (one 
to three percent error In energy) the Relssner-Hxndlln theory can 
replace the equations of the 3-dlmenslonal theory of elasticity at 
least for homogeneous Isotropic plates. The case of anisotropic plates 
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had not been Investigated. The studies have shown that solutions obtained 
from the Relssner-Mlndlln theory are not sensitive to the shear factor. 

5.2 The shear factor 

The most complete engineering analysis of the shear factor to date 
has been presented by Cowper [22]. The shear factors proposed by various 
authors for rectangular sections range between 0.667 to 0.870. The most 
commonly used value Is 0.833 (5/6). The differences arise from various 
corrections for the obviously crude approximation of the shear distri- 
bution Inherent :.ii the Relssner-Mlndlln theory. In his analysis Cowper 
Integrated the equations of the three-dimensional theory of elasticity 
to obtain a formula for the shear factor which accounts for Poisson's 
ratio: 


10(l+v) 

12+llv 


(5.2.1) 


A recent mathematical analysis of the shear factor problem was 
obtained by Vogellus [23]. This analysis Is based on asymptotic expan- 
sion (with respect to plate thickness) of the equations of elasticity 
and accounts for transverse load variations with wave length of the 
order of the plate thickness. Neglecting the effect of rapid trans- 
verse load variations, Vogellus* formula Is: 


. 20(H-v) 

V 24+15V 


(5.2.2) 


We note that for v”0 both formulas yield the same shear factor. For 

V ■ 0.3 K = 0.850 whereas K = 
c V 


0.912 
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5.3 Example; Infinite plate strip, imposed shear displacement 

Let us consider an Infinitely long plate strip of width L and thick- 
ness d subjected to imposed shear displacement, as shown in Fig. S.3.1. 

This problem was chosen for the purpose of studying the effect of shear 
factor on the error of the plate bending solution in relation to the 
plane strain solution. 

The problem can be viewed as a plane strain problem and also as a 
plate bending problem. The finite element meshes and principal boundary 
conditions are shown in Fig. 5.3.1. The plane strain solution was obtained 
with 6 finite elements, quadratlcally graded toward the support, and 
p * 8. Advantage was taken of the two axes of antisymmetry: only one 
quarter of the domain was modelled. The number of degrees of freedom 
for the quarter domain was 399. 

The strain energy of the elasticity solution was estimated by extra- 
polation on the basis of eq. 3.1.3. The results for various L/d ratios 
are shown in Table 5.3.1. 


Table 5.3.1 

Strain energy per unit length of the plane strain solution (U„) . 

E 

6 elements, p = 8, E » 1.0, v = 0.3, L » 2. (One half of the plate only). 


L/d 

2 

5 

10 

20 

40 


0.7648 X 10 
0.7846 X lO”^ 
1.0702 X 10"^ 
1.3668 X 10"^ 
1.7163 X 10"^ 
1.3744 X 10"^ 


200 
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The plate bending model consisted of two finite elements covering half 
of the plate domain, as shown In Fig. S.3.1.b. The polynomial order was 
3, the number of degrees of freedom was 26 for half of the plate domain. 
(This compared with 399 degrees of freedom used for the quarter domain In 
the plane strain solution). The plate bending solution did not change as 
p was Increased to 4, 5 etc. 

The ratios of strain energies computed from the plate bending model 
based on the Relssner^Undlln theory (U„) and from the plane strain 

O 

model (U„) are plotted against the length to thickness ratios for three 
different shear factors In Fig. 5.3.2. It Is seen that for large L/d 
ratios Ug/Ug Is very close to one (the rexiixwc error Is nearly zero) 
for all shear factors. For small L/d ratios the shear factor proposed 
by Vogellus [23] gives the best results. Even when L/d Is only 2, 
the relative error corresponding to Vogellus* shear factor Is less 
than one percent. A large number of studies were performed on the 
effectiveness of the shear factor In more complicated loading situa- 
tions as well; for example: oscillatory loading with the wave length 
of oscillations being equal to the plate thickness. In all cases the 
results were similar to the results presented here. 

The relative errors In maximum deflection, moment and shear force 
were found to be similar In magnitude to the relative error In strain 
energy. 

This example shows that solutions obtained via the Relssner-Mlndlln 
plate theory and the p-verslon of the finite element method are efficient 
and sufficiently accurate for engineering purposes for a very wide range 
of length to thickness ratios. Vogellus* shear factor Is preferable 
when the length to thickness ratio Is less than 20. 
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•-/d 



Fig. 5.3.2 

Infinite plate strip. Imposed shear displacement. 
Ratio of strain energies computed from the plate bending 
model based on the Reissner-Mindlin theory (U ) 
and from the theory of elasticity (Ug) vs. length to 
thickness ratios. 


Relative Error: — — - (percent) 
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5.4 Simply supported square plate under uniform load 

We shall compare the convergence characteristics of the h and p-verslons 
of the finite element method on the basis of the simply supported square 
plate under uniform load, the solution of which Is very smooth. The 
bases of comparison are: L/d ratios and shear factors. 

The rectangular plate domain and the Initial element mesh are shown 
In Fig. 5.4.1 



Fig. 5.4.1 

Simply supported square plate 
Initial finite element mesh 

Figures 5.4.2, 5.4.3 and 5.4.4 illustrate the effect of the length 
to thickness ratio on the rates of convergence. It is seen that as the 
length to thickness ratio Increases from 50 to 1000, the rate of conver- 
gence of the p-verslon based on uniform mesh refinement gradually 




NUMBER OF DEGREES 
OF FREEDOM N 



Fig. 5.4.2 


Uniformly loaded, simply support square 
plate. Length to thickness ratio : 50; 
shear factor: 5/6, Poisson's ratio : 0.3. 


Relative Error (percent) 
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number OF DEGREES 
OF FREEDOM N 



In N 


Fig. 5.4.3 

Uniformly loaded, simply supported square plate. 
Length to thickness ratio : 200 
Shear factor ; 5/6, Poisson's ratio : 0.3 


Relative Error 
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number OF DEGREES 
OF FREEDOM N 



In N 

Fig. 5.4.4 

Uniformly loaded simply supported square plate 
Length to thickness ratio : 1000 
Shear factor : 5/6 Poisson's ratio : 0.3 


Relative Error (percent) 
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number OF DEGREES 
OF FREEDOM N 



In N 

Fig. 5.4.5 

Uniformly loaded, simply supported square 
plate. Length to thickness ratio : 50 
Shear factor : 1000 Poisson's ratio : 0.3 


Relative Error (percent) 
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dlmlnlshes in the preasymptotlc range for p * 2. This Is similar to the 
case of nearly Incompressible solids discussed In conjunction with the 
parabollcally loaded square panel In Section 3.2. For larger p values 
the error la already only about one percent, which Is due to the fact 
that this example problem has a very smooth solution. It Is seen that 
for p " 2 the preasymptotlc rate of convergence Is not significantly 
affected by length to thickness ratios In either the h or the p-verslon. 

The effect of shear factor Is shown In Figures 5.4.2 and 5.4.5. 

In both cases the length to thickness ratio Is 50, however In Fig. 5.4.2 
the shear factor is the commonly used value of 5/6 whereas In Fig. 5.4.5 
the shear factor Is 1000, a value often used to "suppress" shear deforma- 
tion and obtain a solution which Is close to the Klrchhoff plate bending 
solution. The high value of the shear factor has the same effect as 
the high length to thickness ratio: the preasymptotlc rate of conver- 

gence of the h-verslon slows considerably for p * 2. The reason for 
the similarity Is that In both cases the shear strain energy Is multi- 
plied by a penalty term: In one case the penalty term Is the length 
to thickness ratio squared. In the other case It Is the shear factor 
times the length to thickness ratio squared. 

5.5 The rhombic plate problem 

In the preceding section we examined a problem the solution of 
which was very smooth and employed uniform mesh refinement. We shall 
now turn our attention to a much more difficult problem which has 
received a great deal of attention [24,25,26]. In this case we shall 
employ a non-quasiuniform mesh refinement which is "natural" for this 
problem. The problem Itself Is a simply supported, uniformly loaded 
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rhomblc plate, the acute angles of which are 30 degrees. The plate 
and the mesh refinement for one quarter of the plate are shown In 
Fig. 5.5.1. 



Fig. 5.5.1 

The rhombic plate problem 
Non-quaslunlform mesh refinement 

A very strong singularity occurs at the obtuse vertex when 
Klrchhoff *8 theory Is applied to this problem. A detailed discussion 
Is available In [26]. As a result, finite element models employing 
exactly and minimally conforming finite elements fall to converge 
within an acceptable range of number of degrees of freedom. For 
example, Sander reported approximately 17 percent error In the displace- 
ment of the centroid when a conforming displacement model was employed 
with 1000 degrees of freedom [25]. As yet unpublished numerical 


-50- 


number OF DEGREES 
OF FREEDOM N 



In N 

Fig. 5.5.2 

Rhombic Plate Problem 

Relative error In strain energy vs. Number 
of degrees of freedom. Non-quaslunlform 
mesh refinement. Poisson's ratio : 0.3» 
shear factor : 5/6 
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studles performed at Washington University also indicated that fifth 
order exactly and minimally conforming finite elements converge 
very slowly even when rational functions are Introduced. (Detailed 
discussion of rational functions is available in (71). These 
observations underline the importance of employing the Relssner-Mindlin 
theory for modeling structural plates. The severe pollution in 
Kirchhoff plate models caused by obtuse corners is not acceptable 
in practice. 

In this example the finite element mesh is strongly graded toward 

the obtuse corner. As the number of finite elements takes on the 

values n«l, 3, 5, 7... the hypotenuse of the smallest element is 
2 IT — 

10(l-cos ‘^) 2 . In other words, the size of the smallest element 

2 ^ 

decreases in geometric progression with common factor (1 - cos^ 

^ 0.07. Such non-quasiuniform refinement was shown to be optimal when 
coupled with non-quasluniforu p-dlstrlbutlon [9,27]. In this case the 
p-dlstrlbutlon is uniform. The result are shown in Fig. 5.5.2. The 
following observations can be made: 

1. If only one element is used and the relative error in energy 
is to be reduced to under three percent then p must be higher than 8 
whir.h is the hlghast value currently permitted by COMET, (p refers 
to each of the three fields: one transverse displacement and two 
rotations) . 

2. A reasonable grading of elements, such as the 3-element 
mesh, appears to be the best policy: with p ■ 5 only, the relative 
error in energy is under one percent. 
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3. The 5 and 7 element meshes do not provide significantly 
greater accuracy than the 3-element mesh. The reason for this Is 
that the error Is controlled by the largest element, the size of 
which Is not changed In this non-quaslunlform mesh refinement. 

Only an adaptive mesh refinement and/or p-dlstrlbutlon scheme Is 
capable of eliminating the wasteful assignment of degrees of free- 
dom which almost Invariably occurs when the relative contribution 
of elements to the total error of approximation Is unknown. 

This example shows that even In the presence of very severe 
corner singularities It Is possible to obtain approximations, 
accurate to within one percent, with very few finite elements and 
low pol 3 momial order (p = 5) using the Relssner-Mindlin theory of 
plates. The solutions presented here for the 3-element mesh and 
progressively higher p-values are by far the most efficient of 


the solutions published to date for this example problem 
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6. LOCAL ERROR INDICATORS 

Not counting Input errors, there are four sources of error in finite 
element analysis: (1) low polynomial order; (2) singularities due to 

reentrant corners or sudden changes from one type of support condition 
to another along a geometrically smooth boundary; (3) singularities 
due to applied load (e.g. concentrated force) and (4) roundoff errors. 

Low polynomial order can be the source of two kinds of error: it 
may limit the rate of convergence as predicted by eg. (3.1.2) when 
a > p. This is not very important in practical applications because 
a is generally smaller than p. On the other hand, low p values can 
cause non-robust behavior, which may be very difficult to detect In 
practical computations. 

Singularities due to reentrant corners are the most common and 
the least avoidable sources of error. In most practical problems 
there are several corners. In fact, the more complicated the geometry, 
the more frequently will corners of the type discussed in Section 3.1 
occur. Sudden changes in material or geometric properties also belong 
in this category. 

Imposition of severe singularities due to applied loads is generally 
avoidable. The use of concentrated forces in elasticity is a convenience 
but should be avoided when a finite element model based on the displace- 
ment formulation is used because the strain energy associated with 
concentrated forces is unbounded. 

Roundoff errors are less of a problem in the p-version than in 
the h-verslon. The roundoff error depends on the choice of basis 
functions and is, therefore, controllable to a certain extent in the 
development stage of finite element software systems. Full analysis 
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of the roundoff characteristics of the p-version is not yet available. 
The results of a preliminary study were presented in [10]. 

From the foregoing discussion it is evident that in the case of 
the p-version the main sources of error are comer singularities. The 
asymptotic rate of convergence is governed by the most severe of the 
corner singularities. The error depends on the mesh, the polynomial 
order and the loading as well. The contribution of elements to the 
total error of approximation generally varies from element to element. 

In order to minimize the number of degrees of freedom needed for 
achieving a given level of precision, it is necessary to vary the 
p-dlstrlbution over the finite element mesh so as to make the contri- 
bution of each element to the total error approximately the same. In 
the following we shall be concerned with this problem. 

Reliable local indicators for the contribution of individual finite 
elements to the total error in energy are available for the h-version. 
These indicators are based on the residuums and jump discontinuities 
and have the property that the sum of the Indicators estimate the 
total error in energy. 

Previous numerical studies have shown that estimators can be con- 
structed for the p-version as well which have approximately the same 
rate of convergence as the total error in energy but the question of 
whether the same or similar estimators would serve as local Indicators 
for the relative contribution of individual finite elements to the 
total error in energy has not been examined. 

We have investigated a variant of the estimator given in [28], the 
functional form of which was proposed by Babuska [29]. First we define 


our local indicators and demonstrate that the sum of the local indicators 
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Indeed has the same rate of convergence as the error In energy when uniform 
p-dlstrlbutlon Is used on a fixed finite element mesh. We then demonstrate 
through example problems that the indicators also serve to estimate the 
relative contribution of individual finite elements to the total error 
in energy and provide for establishing optimal or nearly optimal p-dlstributlon 
for a given mesh. Finally, alternative policies for achieving optimal 
p-dlstributlon in practical computations are discussed. 

6.1 Definition of local error indicators 

Local error indicators are computed from the reslduums and jump 
discontinuities of the finite element solution. The reslduums r^ are 
obtained from substituting the finite element solution into the equilibrium 
equation. In elasticity: 

r. + (6.1.1) 

in which r^ is the i^^ component of the unbalanced body force (residuum) ; 

X and G are the Lame parameters, TT is the i^^ component of the displace- 
ment vector computed by the finite element method; is the Imposed 
body force. 

The jump discontinuities are the differences between the surface 
tractions computed from the displacement functions belonging to 
neighboring elements along the common finite element boundary. We shall 
denote jump discontinuities by At^^^ in which 1 is the edge number, m is 


the usual component index 
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Our goal Is to construct an a posteriori error indicator which measures 
the relative contribution of each finite element to the total error of 
approximation. The error indicator for the element consists of 
two parts; one based on residuums, the other based on jump discontinuities, 
and has the following form: 




+ C T, 


( 6 . 1 . 2 ) 


in which is the *'R indicator"; C is a constant; is the "T indicator" 
for the k^^ element. We shall now briefly outline the rationale for 
constructing Rj^ and Tj^. 


6.1.1 The R indicator 

The residuum computed from the finite element solution is propor- 
tional to the error in the second derivatives: Suppose that the exact 

solution u. is known. Then: 

1 

° ^ '"i.jj “j,ji h (6.1.3) 

and, from equations 6.1.3 and 6.1.1: 


We seek to approximate the error in strain energy which is the integral 
over the finite element domain of a linear combination of the squares of 
the error in the first derivatives. In view of this fact, the R estima- 


tor is constructed as follows: 
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5 


ij ri rj dA + 


(6.1.5) 


in which 


E 

Pk 


r 


i 


6 


ij 


W,. 


Is the area of the k'^'^ element; 
is the modulus of elasticity; 

Is the polynomial order of the displacement approximation 
over the k^^ element; 

Is the unbalanced body force vector component, defined 

In eq. 6.1.1 

Is the Kronecker delta; 

Is a weight function, defined as the square root of the 

first Internal mode; l.e. the square root of a cubic 

polynomial that vanishes along the boundaries of the k^^ 

element. W, Is normalized such that: 
k 


IT I 


dA = 1 


( 6 . 1 . 6 ) 


Rj^ and Rj^ are correction terms which are zero when the p-dlstrlbutlon 
over the k^^ element and Its Immediate neighbors Is uniform. When this 
condition Is not met then the k^^ element has higher and/or lower order 
neighbors than Pj^. The lower order neighbors downgrade the k^^ element 
along their common boundaries. Similarly, the k^^ element downgrades 
Its higher order neighbors. An Indicator of the relative contribution 
of the k^^ element to the total error In energy must In some way account 
fur these effects. R^ and ^ were devised for that purpose on the basis 
of extensive numerical experiments. We shall now define these correction 


terms: 


-58- 


Let be the set of element numbers which have common edge with 
element k. Also, let 


= min {p^} 

iel,. 


(6.1.7) 


\ Pi < Pk . i I 


If pj^ > pj^ then Rj^ = 0. Otherwise: 


4-2V'-L-”^ f 

2 - 2 " 2 El 

±e\ J J \ 


6 r r W, dA 
mn m n k 


(6.1.8) 


Similarly, let pj^ = max {p^} and 


= {i Pi > Pjt ’ ^ ^ ®k - ^k ^ 


Otherwise: 


V' _L. f X u 

■ y. 2 2^ 2 - 2 " 2 El *^001 ’^m *^n 

^ Pi _ P4 Pi Pi J 

iel jel J ^ ^ A 

^ (6.1.9) 


is the correction term that accounts for the downgrading of the k^ 
element by its neighbors; ^ accounts for the k*"^ element downgrading its 
neighbors. 


Further details are available in reference 30. 
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6.1.2 The T Indicator 

The T-lndicator accounts for jump discontinuities at interelement 
boundaries and the error between computed and specified tractions at 
external boundaries: 

Let be the set of edges of element k for which surface tractions 
* 

are specified and let ^k' 



in which 


til 

is the length of the i edge of element k; 

'’ik ■ Pfc + Pi ^ P -k ! 

= 2 Pj^ when i e and no displacement vector component 

is prescribed along edge i; 

Plk ® Pj^ when i e and one displacement vector component is 

prescribed along edge 1; 

At is the difference between the surface tractions 

m 

calculated from the finite element solutions for 
element k and element i along their common edges 
when i e I, . When i e I, then At ' is the difference 
between the surface tractions calculated for element k 
and the prescribed surface tractions along edge 1; 
is the square root of the first Internal edge mode 
for edge 1, normalized such that 
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-j- f d* . 1 

1 

Finally, we define the error Indicator for the element as: 

■ V * ‘ \ (6.1.U) 

In which c = - with summation over all elements. 

6.2 Example: Parabollcally loaded square panel 

This problem Is the same as the problem discussed In Section 3.2, 
however the finite element mesh Is different (see Fig. 6.2.1). 



Fig. 6.2.1 

Parabollcally loaded square panel 
Finite element mesh and element numbers 
Poisson's ratio : 0.3 
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We first demonstrate that and (summation over all elements) 

have the same rate of convergence as the total error In energy when 

uniform p-dlstrlbutlon Is used on a fixed finite element mesh. In 

Fig. 6.2.2 we have plotted the error In strain energy (l.e. the error 

2 

In energy norm squared | | e | | ^) ; and on log scale against 

N also on log scale. It Is seen that the slopes of the three curves 
are very nearly Indentlcal, especially for large N, Indicating that the 
asymptotic rates of convergence are also nearly Identical. 

Next we demonstrate that Is a reliable Indicator of the relative 
contribution of element k to the total error In strain energy: Starting 

from uniform p « 2, we Incremented p, one element at a time by one, for 
all elements, and computed the corresponding errors In energy. We then 
selected that p-dlstrlbutlon for which the error was the smallest and 
using this p-dlstrlbutlon as the base we again Incremented p by one, 
one element at a time, for each element and selected that p-dlstrlbutlon 
for the next base for which the error was the smallest. We continued this 
process of exhaustive evaluation until p had to be Incremented beyond 8 
which Is the limiting value for our computer program, COMET-X. In this 
way we obtained an optimal sequence of p-dlstrlbutlons. Next, starting 
again from uniform p » 2, we computed ^or all elements and Incremented 
p by one In that element only for which was maximum. Proceeding In 
this way we obtained an Indicated sequence of p-dlstrlbutlons. Both 
the optimal and Indicated sequences are shown In Table 6.2.1. The two 
sequences are not Identical only In the third and fourth cycles. For 
these cycles the results of the exhaustive evaluation are tabulated. 

Under "Remarks" the Indicated and optimal sequences are noted. The 



Logdlellg, IR,j, IT,,) 


P 



Log {'/(yj) 

Fig. 6.2.2 

Parabolically loaded square panel 
Error in strain energy and sums of the R and T 
indicators over all finite elements 
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Table 6.2.1 

Optimal and Indicated Sequences of P-Distrlbution for 
the Parabolically Loaded Square Panel 


p 

N 

u-u 

100 ^ 

P 

Remarks 

2 2 2 2 

20 

1.70 

0.53 


3 2 2 2 

25 

6.68 E-1 

0.48 


4 2 2 2 

32 

5.38 E-1 

0.90 

Optlioal 

3 3 2 2 

30 

5.66 E-1 

0.80 


3 2 3 2 

29 

6.51 E-1 

0.79 

Indicated 

3 2 2 3 

29 

6.68 E-1 

0.80 


4 2 3 2 

36 

5.16 E-1 

1.12 

Indicated 

3 3 3 2 

36 

1.13 E-1 

0.64 

Optimal 

3 2 4 2 

33 

6.50 E-1 

1.02 


3 2 3 3 

35 

6.51 E-1 

1.15 


mm 

43 

6.15 E-2 

0.74 


mm 

50 

1.63 E-2 

0.61 


4 4 4 2 

58 

4.34 E-3 

0.51 



67 

2.94 E-3 

0.59 



Note: The polynomial orders are listed in the order of element 

numbering shown in Fig. 6.2.1. 
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Table 6.2.2 

Sequence of P-Dlstrlbutlons Obtained on the Basis of Elemental Errors 
In Strain Energy for the Parabollcally Loaded Square Panel 


p 

N 

100 

P 

2 2 2 2 

20 

1.70 

0.53 

3 2 2 2 

25 

6.68 E-1 

0 48 

4 2 2 2 

32 

5.38 E-1 

0.90 

5 2 2 2 

41 

5.29 E-1 

1.47 

6 2 2 2 

52 

5.29 E-1 

2.36 

7 2 2 2 

65 

5.28 E-1 

3.68 

8 2 2 2 

80 

5.28 E-1 

5.58 


See footnote to Table 6.2.1 






-65- 


number of degrees of freedom, the relative error in energy (percent) and a 
cost ratio p, defined as the approximate cost associated with the non- 
uniform p-dlstrlbutlon divided by the corresponding cost for uniform 
p-dlstrlbutlon for the same relative error In energy (assuming asymptotic 
behavior) are also shown In Table 6.2.1. 

Occasional differences between the Indicated and optimal sequences 
can be explained by pointing out that the effect of downgrading one 
element by Its neighbors cannot be estimated with precision. It Is 
Important only that the Indicators be stable In the sense that If 
deviations occur from the optimal path then, following the Indicators, 

It should be possible to return to the optimal path. We note also that 
small deviations are unimportant In practical computations where p would 
be Incremented concurrently for groups of elements for which the Indica- 
tors are within a certain tolerance range. 

Finally, we make the following observation: The Indicators are not 

to be confused with Indications of error In strain energy for Individual 
finite elements. In fact, the error In strain energy for Individual 
finite elements Is a very poor Indicator. To Illustrate this, starting 
from uniform p ■ 2 we incremented p by one at a time for that finite 
element only for which the error In strain energy was the largest. 

The resulting sequence of p-dlstribution Is shovm in Table 6.2.2. The 
sequence Is rapidly deviating from the optimal one which can be Judged 
from the progressively Increasing value of p. 

6.3 Example: Double-edge cracked square panel 


Our next example Is the double-edge cracked square panel, shown In 
Fig. 6.3.1. The problem Is essentially the same as the problem discussed 
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<r 


X 


Fig. 6.3.1 

Double-edge cracked square panel 
Finite element mesh and element numbers 
Poisson's ratio : 0.3 

under section 4.3, but the finite element mesh is different. As in the 
preceding example, we obtained an optimal sequence of p-distributions 
by exhaustive evaluation, however in this case we started from p 3 
(uniformly). Also, we obtained an indicated sequence of p-distributions 
using the indicator The results are given in Table 6.3.1. Remarkably, 

the Indicated sequence deviates from the optimal one only in the seventh 
and eighth cycl s. The relative error for the optimal and indicated 
sequences is shown in Table 6.2.1 and plotted in Fig. 6.3.2. 
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Table 6.3.1 

Optimal and Indicated Sequences of P-Distribution for the 
Edge Cracked Square Panel 










































U(|i)*U(UrE) 



Fig. 6.3.2 


Double-edge cracked panel 
Relative error vs. number of degrees of 
freedom for the optimal and indicated 
sequences of p-distribution 


Relative Error (Percent) 
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Thls example shows that the error indicator is a reliable guide even 
when a strong singularity is present. It is also seen that the cost savings 
achievable with optimal p-distrlbutlon, as compared with uniform p-dlstribu- 
tion, are between 40 and 50 percent (the value of 1-p expressed as percen- 
tage) . These savings are strongly mesh dependent, however. For finer 
meshes the savings would be greater; for coarser meshes the savings would 
be smaller. 

6.4 Example: Lap joint problem 

In this problem we have a fairly complicated geometry and 10 geometric 
singularities of various Intensity (Fig. 6.4.1). The plates are Joined 



Fig. 6.4.1 
Welded lap-joint 

Finite element mesh and element numbering 
Poisson's ratio : 0.3 
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Table 6.4.1 

Optimal and Indicated Sequences of P-Dlstributlon 
for the Lap-Joint Problem 


p 

— 

N 

100 

P 

Remarks 

2 2 2 2 2 2 2 

72 

29.539 

21.22 


3 2 2 2 2 2 2 

80 

25.935 

20.20 


3 3 2 2 2 2 2 

92 

4.418 

0.77 


4 3 2 2 2 2 2 

104 

3.688 

0.69 


4 3 3 2 2 2 2 

— 

3.244 

0.62 

Optimal 

4 4 2 2 2 2 2 


3. 435 

0.80 

Indicated 

4 3 3 3 2 2 2 


2.896 

0.60 

Optimal 

5 4 2 2 2 2 2 


3.128 

0.85 

Indicated 

5 3 3 3 2 2 2 

140 

2.579 

0.61 

Optimal 

5 4 2 3 2 2 2 

144 

2.832 

0.78 

Indicated 

5 4 3 3 2 2 2 


2.283 

0.59 

Optimal 

5 5 2 3 2 2 2 

■91 

2.536 

0.81 

Indicated 

5 5 3 3 2 2 2 

176 

2.050 

0.61 


5 5 3 3 3 2 2 

188 

1.850 

0.57 

Optimal 

5 6 3 3 2 2 2 

196 

1.955 

0.69 

Indicated 

6 5 3 3 3 2 2 

208 

1.691 

0.58 

Optimal 

5 6 3 3 2 2 3 

204 

1.839 

0.66 

Indicated 

6 5 4 3 3 2 2 

220 

1.554 

0.54 

Optimal 

6 6 3 3 2 2 3 

228 

1.638 

0.65 

Indicated 

6 6 4 3 3 2 2 

244 

1.416 

0.56 

Optimal 

6 7 3 3 2 2 3 

252 

1.554 

0.72 

Indicated 

6 6 4 3 3 2 3 

256 

1.321 

0.54 

Optimal 

7 7 3 3 2 2 3 

280 

1.469 

0.79 

Indicated 

6 6 4 3 3 3 3 

272 

1.162 

0.47 

Optimal 

7 7 3 4 2 2 3 

292 

1.448 

0.84 

Indicated 

7 6 4 3 3 3 3 

296 

1.088 

0.49 

Optimal 

7 7 3 4 3 2 3 

308 

1.268 

0.72 

Indicated 

7 7 4 3 3 3 3 

■9 


0.50 

Optimal 

7 7 4 4 3 2 3 

B9 


0.62 

Indicated 

7 7 4 3 3 3 4 

336 

0.941 

0.47 

Optimal 

7 7 4 4 3 3 3 

340 

0.951 

0.49 

Indicated 

7 7 4 4 3 3 4 

352 

0.898 

0.47 
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number OF DEGREES 
OF FREEDOM (N) 



In N 

Fig. 6.4.2 
Welded lap-joint 

Relative error vs. number of degrees of freedom 
for the optimal and indicated sequences 
of p-distribution 


i 

1 
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only by Che first welds (element 3) and are free at the base of element 7. 
Loading Is by imposed displacement, as shown in Fig. 6.4.1. 

The sequences of p-dlsCributlon (arranged in the order of element 
numbering) obtained from exhaustive evaluation and from evaluations based 
on Che error indicator are shown in the first column of Table 6.4.1. The 
relative error is plotted against the number of degrees of freedom for 
both sequences in Fig. 6.4.2. 

The results show that the indicated sequence does not diverge from 
the optimal sequence. As previously noted, local deviations are attribu- 
table to the fact that the effect of dotmgrading of one element by its 
neighbors cannot be estimated with precision: depending on Che approxi- 
mated function, Che suppressed degrees of freedom may be very important 
or not Important at all from the point of view of contribution to Che 
total error of approximation. Numerical experiments have indicated 
that deviations are less apt to occur when p is Increased in more than 
one element at a time: for example p could be Increased for all elements 

simultaneously for which the error indicator is greater than the average 
value of all error indicators. 
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7. STRUCTURAL SYNTHESIS 

In the preceding sections we focused on the reliability and efficiency 
of the analysis part of the programming system for structural synthesis. 

We have shown that the desired level of precision closely determines 
the required number of degrees of freedom, hence the cost of analysis, 
which Is the dominant part of the structural synthesis process from 
the point of view of resource requirements. 

The phrase: "desired level of precision" means two things: the 

numerical value of the acceptable relative error and definition of the 
norm In which the error Is to be measured. In our analysis the error 
was measured In energy (the square of the energy norm In use In mathe- 
matical works). In a certain sense this Is the logical choice because 
the finite element method based on the displacement formulations 
actually minimizes the error In energy norm. Other norms, often 
regarded as more convenient for engineers, do not always exist: For 

example, In the case of the lap-joint problem, discussed In Section 6.3, 
there are 10 geometric singularities, 8 of which result In Infinitely 
high stresses. In this case It would not be meaningful to state the 
desired level of precision In terms of stresses: The error Is Infinitely 

large at the singular points. The stresses are approximated In the least 
squares sense, rather than polntwlse, In the displacement formulation of 
the finite element method. 

In structural synthesis the constraints may be stated In a variety 
of different norms: stress resultants (moments, axial and shear forces) 
at specific points, displacements, vibration frequencies, margins of 
safety against buckling, etc. In general, the corresponding acceptable 
tolerance levels can be given on the basis of engineering considerations. 
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In order to minimize the cost of analysis It Is necessary to select 
the number and distribution of the degrees of freedom through proper 
choices of the mesh and polynomial orders so that the acceptable 
tolerances are neither exceeded nor satisfied with excessive margins. 

Ideally, one would need reliable local a posteriori error estimators 
In the various norms to achieve this goal. At present such estimators 
are not available. Some preliminary work, has been completed under 
the present project for the p-version in energy norm only, which was 
reported In Section 6, but much fundamental work remains. The state 
of the art Is more advanced for the h-verslon [28] but the computer 
Implementation Is more difficult. The benefits to be gained from 
such estimators are twofold: Increased confidence In the computed 

results and reduced cost. The cost reduction will be problem depen- 
dent, however It can be estimated to be about 50 percent for most 
practical problems, based on what Is already available through appli- 
cation of the p-verslon of the finite element method. Reliable local 
a posteriori estimators In other norms useful In engineering analysis 
can be developed In two to five years. Full utilization of adaptivity 
In structural synthesis requires the analyzer to communicate to the 
solver the required tolerance levels at the beginning of each Iterative 
cycle. 

Even In the absence of fully adaptive selection of the proper 
p-dlstrlbutlon, the p-version provides for very substantial cost reduc- 
tions In the analysis part of the structural synthesis process. The 
basis for cost reduction Is the faster rate of convergence of the 
p-verslon which can be realized almost Invariably In engineering analysis. 
(Wave propagation problems In which the wave front Is not at finite 
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elemenc boundaries are exceptions: the asymptotic rates of convergence 

of the h and p-verslons In such cases are Identical [9,10]). 

The rationale for utilizing Indirect error estimation techniques 
In conjunction with the p-verslon of the finite element method and 
the recommended procedure is as follows: It Is known from theoretical 

analyses that the rate of convergence (l.e. the relative error vs. the 
number of degrees of freedom relationship) Is governed by geometric 
and loading singularities. This was discussed In Section 3 of this 
report. Additional details are available In references 9 and 10. The 
various types of singularities can be readily classified on the basis 
of the smoothness parameter of the approximated function as It was done 
by Williams [12] for plane elastic and plate problems. A series of 
benchmark problems, characterized by various smoothness parameters, 
can now be used for establishing the relationship between relative 
error and polynomial order. In a given problem, where several 
different singularities occur, the conservative approach Is recommended 
the polynomial order should be selected on the basis of the strongest 
singularity. For example, referring to figure 4.3.2 we see that 
p»5 corresponds to a 2 percent relative error for the edge cracked 
panel problem using a coarse mesh. In fact. If only a 3-element mesh 
were used the relative error would still be about two percent for 
P-5 [31]. 

The smoothness parameter, which measures the strength of the slngu 
larlty, (the value of a In eq. 3.1.1) Is 0.5 In this case. Similarly, 
the strongest singularity in the welded lap-joint problem has the 
same smoothness parameter: a > 0.5. Thus, If we wish to have a 
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solutlon which Is accurate to approximately 2 percent In energy for the 
welded lap-joint problem, we should employ p ■ 5 uniformly. 

7.1 Example; Simply supported rectengular plate with moveable supports . 

In this example we compare the computational costs of optimizing 
a simple problem via the h and p-versions of the finite element method 
requiring similar levels of precision for both versions and employing 
existing technology. 

The problem definition is given in Fig. 7.1.1. The plate is simply 
supported, the moveable support is also a simple support. Because of 
symmetry, only one quarter of the plate needs to be modeled. The 
boundary conditions and the four-element mesh are shown in fig. (b), 
the 6 and 96-element meshes are shown in figures (c) and (d) . In order 
to perserve uniformity, in the subsequent discussions the plate thick- 
ness is 1.0 in; Poisson's ratio is 0.2; the modulus of elasticity is 
30,000 ksl, the shear factor is 3/6; the lateral load is uniform 
pressure of 5.0 ksi and the position of the intermediate support (x) 
is 6.0 Inches, unless otherwise noted. 

We note that there are two singularities, one at point A the other 
at point B in Fig. 7.1.1(a). It is possible to compute the strength 
of these singularities by analytic methods, but not necessary. It is 
sufficient to analyze a model problem with the same singularities, to 
establish the minimum p-level and minimum mesh refinement needed for 
the finite element solution. 

Because this is a simple problem, we use the problem itself to 


find the minimal p-level and mesh refinement. This is done in 
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Moveabie 

SuDDort 





(a) 






(d) 

Fig. 7.1.1 


Simply supported rectangular plate with 
uoveable internal supports, (a) Dimensions; 
(b) 4-element mesh and principal boundary 
conditions; (c) 6-element mesh; 

(d) 96-element mesh. 



Strain Energy (U 
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N 



1.0 1.5 2.0 2.5 3.0 

Log N 


Fig. 7.1.2 

Sliaply supported rectangular plate 
with moveable supports. Four 
Element Mesh. Strain energy vs. N. 
Support location : x * 6.0 
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N 

20 100 1,000 10,000 



Log N 


Fig. 7.1.3 

Simply supported rectangular plate 
with moveable supports. Relative 
error in strain energy vs, N 
Support location : x ■ 6.0 


elative Error (Percent 
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two steps: First we compute the approximate value of the true strain 

energy, using the 4-element mesh In conjunction with the p-verslcn of 
the finite element method. The results are shown In Fig. 7.1.2. As 
p Is Increased, the value of the strain energy approaches 7.4054 In-k 
for the quarter plate. In order to Illustrate the fact that this 
result Is not sensitive to the mesh, we did a similar analysis for 
the 6-element mesh and obtained the same result. Using 7.4054 In-k. 
as the estimated true strain energy, we next computed the two exten- 
sions (by the h and p versions) of an Initial solution In the relative 
error vs. number of degrees of freedom diagram, with logarithmic scales 
(Fig. 7.1.3). In this case our Initial solution was the 6-element 
mesh and p = 2. For the h-version p was 2 for all three fields and 
the mesh was refined quasi-uniformly, as shown in Fig. 7.1.1(d). The 
refinement allows for varying the position of the support without 
introducing excessive distortion In the elements. Of course, any 
h-verslon finite element code could have been used for the h-extenslon. 
In the Interest of preserving uniformity we employed COMET for both the 
h and p-versions. 

It Is seen from Fig. 7.1.3 that the 96-element mesh with p 2 
and the 6-element mesh with p » 4 yield about the same relative error 
In strain energy which Is between one and two percent. This level of 
accuracy Is probably more than adequate for most engineering purposes. 
The key computational parameters are listed In Table 7.1.1. 
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Table 7.1.1 

Comparison of solutions by the h and p versions 
on the basis of similar accuracies 


h-version 


3-version 


No. of elements 

96 

6 

Poljmomial order (p) 

2 

4 

Relative error in strain 
energy (percent) 

1.45 

1.19 

No. of degrees of freedom 

600 

156 

Total solution time (CPU sec.) 

251 

41.5 

Solution of equations (CPU sec.) 

61 

14.5 

Computation of stiffness matrices 

125 

15.3 


and load vectors (CPU sec.) 


The computations were performed on a DEC System 2040 computer; 

128K 36-1 it work memory, TOPS-20 Operating System. Because the Indicated 
CPU time varies with the overall workload of the machine, exact dupli- 
cation of the CPU times shown in Table 7.1.1 is not possible. Varia- 
tions can be as high as 10 percent. 

The results are typical in the sense that the total solution time 
for the h version is about 6 times that for the p-version. In our 
structural synthesis studies we varied the position of the support 
(x) and the plate thickness (t). The constraint was that the maximum 
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stress, computed from the moments In the usual way, should not exceed 
20 ksl. For both the h and p versions convergence occurred after 
10 or 11 iterations. Thus, using the h-version, the solver required 
approximately 42 CPU minutes whereas using the p-version the solver 
required approximately 7 CPU minutes. The results of Section 6 
indicate that with the use of optimal p-dlstribution additional savings 
of 50 percent should be possible. 

In extrapolating these results for other problems, the following 
points must be considered: (1) The CPU time vs. N curve does not pass 

through the origin because a certain amount of "overhead" is involved in 
initiating each solution. (2) The cost of computation of elemental 
stiffness matrices and load vectors increases linearly with the 
number of elements and with the square p. (3) The cost of solution of 
the system of simultaneous equations increases approximately with the 
square of the number of degrees of freedom. For very large problems 
this cost overwhelms all other costs, (4) The cost ratio between h- 
and p-version solutions depends on the strongest singularity when 
quasiuniform meshes and p-distributlons are used. In general, the 
cost ratio is more favorable for the p-version when the singularity 
is stronger. In this example the singularity was relatively mild. 

Had the problem contained an obtuse corner as in the case of the rhombic 


plate discussed in Section 3.5, the cost ratio would have been extremely 
high. 



8. CONCLUSIONS AND RECOMMENDATIONS 


Adaptive finite element technology has the potential for reducing 
the cost of structural stress analysis by about an order of magnitude. 

This will strongly Impact structural synthesis where one of the current 
limitations Is the cost of analysis. 

Adaptivity can be based on either the h- or p-verslons of the finite 
element method, or a combination of both. From the practical point of 
view, the goal being to obtain reliable approximations accurate to within 
one and three percent relative error In energy, adaptivity based on the 
p-version Is the most promising approach. In particular, the p-verslon 
offers greater reliability through its robust behavior and greater efficiency 
(due to the fact that progressive mesh refinement Is not required) than the 
h-verslon. Within the accuracy range required in engineering analysis, 
combination of the h> and p~verslons does not offer practical advantages. 

Full utilization of the adaptive process in structural s)mthesis 
requires that local a-posteriori error estimators be developed in those 
norms in which the constraints are specified. At present only preliminary 
results are available for the p version and only in energy norm. The 
state of the art is more advanced in the case of the h-version. We 
recommend that research should be vigorously pursued in this area and 
estimate that local a posteriori error estimators can be developed in 2 
to 5 years. The estimated payoff would be approximately 50 percent 
reduction in cost over presently available adaptive approaches. 

Certain adaptive strategies are already available. We refer to 
the fact that reasonably designed fixed meshes, coupled with uniform 
p-distribution, with p sufficiently high to guarantee the desired level 
of precision, can already result in very substantial reductions in the 
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cost of analysis. The amount of reduction Is problem-dependent, nevertheless 
for practical problems It can be as high as several orders of magnitude, 
when compared with conventional finite element technology. 

The choice of p-level must be based on considerations of the desired 
level of precision. In the absence of suitable a posteriori estimators. 
Indirect error estimation should be used. Indirect error estimation Is 
based on the fact that the parameters controlling the accuracy of solution 
In the p-verslon can be determined before starting the synthesis process. 

An Important advantage of the p-verslon, not previously realized. Is 
Its robustness. Examples were presented which Indicated that the p-verslon 
Is Insensitive to Poisson's ratio In the case of elasticity and to plate 
thickness In the case of the Relssner-Mlndlln theory of plates. This 
ensures proper convergence even when the plate or shell thickness ranges 
between wide limits In the structural synthesis process. 

It has been shown that the Relssner-Hlndlln theory of plates Is 
capable of approximating the three dimensional elasticity solution to 
within one or two percent relative error In energy for all length-to 
thickness ratios and loading conditions likely to occur In practical 
applications of structural synthesis. The same statement cannot be 
made of the Klrchhoff theory of plates. For this reason, use of the 
Relssner-Mlndlln theory Is recommended for practical analyses. We 
remark, however, that this conclusion and recommendation are limited 
to plates manufactured from Isotropic materials. Composite plates 
may present convergence problems unless the shear factor Is carefully 


chosen. 
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